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A precise evaluation of the secondary particle production and propagation in the atmosphere 
is very important for the atmospheric neutrino oscillation studies. The issue is addressed with the 
extension of a previously developed full 3-Dimensional Monte-Carlo simulation of particle generation 
and transport in the atmosphere, to compute the flux of secondary protons, muons and neutrinos. 
Recent balloon borne experiments have performed a set of accurate flux measurements for different 
particle species at different altitudes in the atmosphere, which can be used to test the calculations for 
the atmospheric neutrino production, and constrain the underlying hadronic models. The simulation 
results are reported and compared with the latest flux measurements. It is shown that the level 
of precision reached by these experiments could be used to constrain the nuclear models used in 
the simulation. The implication of these results for the atmospheric neutrino flux calculation are 
discussed. 
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PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

The need of a precise knowledge of the atmospheric 
production of secondary particles induced by the cosmic 
ray (CR) flux of charged particles, became a major stake 
for the scientific community with the recently reported 
evidence for atmospheric neutrino oscillations pointing 
to a non-zero neutrino mass, which evidence was based 
on a comparison of the measured data with secondary 
atmospheric neutrino flux calculations. In this context, 
the rapidly increasing amount and statistical significance 
of the data collected by underground neutrino detectors 
PI , also made precise atmospheric neutrino flux calcu- 
lations highly desirable as a potential tool to test calcu- 
lations and models. Onc-dimcnsional calculations of the 
neutrino flux have been considered for long as a good ap- 
proximation of the flux at earth, until some first attempts 
of three dimension calculations were performed 0, 0,11, 
which rapidly became the new standard approach. A 
survey of the different approaches can be found in 0] . A 
precise calculation of the neutrino flux relies on a precise 
knowledge of the hadronic production cross sections, of 
nucleons and mesons at the top of the decay chains lead- 
ing to neutrinos in the final state. These cross sections 
however are not known in general with a high level of 
accuracy. 

A three-dimensional (3D) simulation describing the 
CR induced cascade in the atmosphere, particle prop- 
agation in the geomagnetic field, and interactions with 
the medium, was developed by the authors for the inter- 
pretation of the AMS01 data. The code was successfully 
used to reproduce theproton, electron-positron, and he- 
lium 3 flux data 0, U 01 measured by AMS and their 
respective dependence on the geomagnetic coordinates. 
The lack of a 3D calculation of the neutrino flux at this 
time led the authors to be the first to investigate this 
important issue either, with some first results reported 



in A complete report on the calculated muon and 

neutrino flux was published in [Toj . 

This paper reports the results of a further investigation 
of the issue, with the calculations improved by several 
respects. One improvement consisted in the use of vari- 
ance reduction techniques (see section Til Fl below) which 
allowed, together with some code optimization, to ap- 
preciably increase the statistics of the simulated events 
sample, and thereby improve the Monte-Carlo precision 
(statistical accuracy) of the calculation. A second im- 
provement was based on the observation that the par- 
ticle production cross sections are the main source of 
uncertainty in the secondary particle flux calculations. 
The reliability of the neutrino flux simulation is usually 
tested on the capacity of the calculations to reproduce the 
atmospheric muon and proton flux measurements. The 
latter can be as well used to directly constrain the sec- 
ondary particle production cross sections, and to reduce 
the simulation uncertainties on the calculated neutrino 
flux. This prospect has been explored in details and it is 
shown (see section lTll C|) that the accuracy of the flux cal- 
culations can be significantly improved by this method. 

The paper is organized as follows. The method and 
models used in the calculations are introduced in section 
ITT1 Section II I II is devoted to the results on the proton, 
muon and neutrino flux in the atmosphere. Summary 
and conclusions are given in section llVl 



II. SIMULATION MODEL 

The calculation proceeds by means of a full 3D- 
simulation Monte Carlo simulation. In this section the 
main features of the pro gram are recalled for convenience. 
See refs HH II El El for other details. 
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A. Cosmic ray flux 

For the primary flux, the 1998 AMS measurement of 
CR proton and helium flux ^| an d the fit from 0] 
for heavier elements (up to iron), are used. The kinetic 
energy range of incident CRs covered in the simulation 
is [0.2, 10000] GeV/n. 

For each periods of the solar cycle, the incident cos- 
mic flux are corrected for the different solar modulation 
effects using the simple force field approximation [lfij . 

Incident Cosmic Rays are generated on a virtual sphere 
at a finite distance in the earth neighborhood. The flux 
at any point inside the volume of this virtual sphere is 
isotropic provided the differential element of the zenith 
angle distribution of the particle direction generated on 
the sphere is proportional to cos^ d(cos9 z ), 6 Z being the 
zenithal angle of the particle. The geomagnetic cut-off is 
applied by backtracing the particle trajectory in the geo- 
magnetic field (using the method described in lll C|) , keep- 
ing in the sample only those particles reaching a back- 
tracing distance of 10 Earth radii. The altitude of the 
virtual sphere has to be chosen outside the atmosphere 
to ensure a backtracing process in a region free of inter- 
action. In the present case it was chosen at 2000 km. 

For cosmic ray nuclei, the superposition model |l6j was 
used. In this approximation a nucleus of mass number A 
and charge number Z is replaced by A — Z neutrons and 
Z protons with the same velocity as the parent nucleus. 
The particles are then processed like real protons and 
neutrons, excepted that the effective charge used for the 
propagation in the geomagnetic field is set to Z/A to keep 
the same rigidity (i.e. the same trajectory) as the parent 
nucleus. 



B. Atmospheric and Geomagnetic model 

The model used to simulate the atmosphere is the 
MSISE-90 model available from which describes the 
temperature and densities in the earth's atmosphere from 
ground to thermospheric altitudes. 

The geomagnetic model used in this version of the sim- 
ulation is IGRF9 0| which is the ninth generation of a 
mathematical model of the Earth's main field and its an- 
nual rate of change (secular variation) adopted by the 
International Association of Geomagnetism and Aeron- 
omy (IAGA). 

C. Particle propagation 

Each particle is propagated in the geomagnetic field 
and interacts with nuclei of the local atmospheric den- 
sity The particles trajectory is built by numerically 
integrating the equation of motion using the fifth or- 
der Runge Kutta method with adaptive step-size control 
from ref |23|. At each step of the propagation, in addi- 
tion to the three coordinates and the three components 
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FIG. 1: Proton invariant triple differential cross section mea- 
sured (open square) in the p + Be — > p + X at 19 Gev/c reac- 
tion \u\ compared with the fitted parametrization as a func- 
tion of the proton kinetic energy and for different production 
angles : 0.7°, 1.1° (scaled by a factor 10" 1 ), 1.7° (xlO -2 ), 
2.3° (xl(T 3 ), 2.9° (xlO -4 ), 3.4° (xlO -5 ) , 4.0° (xlO -6 ). 



of the momentum, the grammage crossed by the particle 
is computed. This allows the integration step size to be 
adapted to the local value of the atmosphere density and 
then insure a smooth penetration in the atmosphere and 
a precise calculation of the interaction probability and of 
the ionization energy loss, and the latter be computed 
for each step along the trajectory. 

Every secondary particles are processed the same way 
as their parent particle, leading to the generation of an at- 
mospheric cascade, more or less extensive, for each event. 
Nuclcons, pions and kaons are produced for each interac- 
tion with their respective triple differential cross sections. 



D. Particle interaction and secondary production 

The cross sections used to produce secondary particles 
(p, n, 7r , K^) are a very important input in the cal- 
culation, and the main source of uncertainty in the at- 
mospheric flux estimation. In the present approach the 
particle production cross sections are obtained from fits 



3 



to the data using an approach based on the Kalinovsky- 
Mokhov-Nikitin (KMN) parametrization of the inclusive 
hadronic cross sections j22f . In the present work, a wide 
set of experimental data have been used to constrain the 
parameters of a modified KMN analytical formula |2lj . 
The set of data included in the fits have been selected to 
cover as far as possible, the kinematics of interest for at- 
mospheric secondary particles production. An interest of 
this method is that together with the best fit parameter 
set for each reaction channel, the errors on the parame- 
ters can also be estimated. 

An illustrative example is shown on figure^ with the 
comparison between the data from ^3] (j> + Be — ► p + X 
at 19 Gev/c) and the fit results (solid line). For all the 
studied samples of each reaction channel, the reduced 
X 2 obtained from fits were found larger than 1 , typically 
between 1 and 5. This could be expected since the sta- 
tistical and systematic errors on the cross section mea- 
surements arc usually approximate or poorly known. To 
account for this uncertainty, a scale factor on the ex- 
perimental error equal to l/\Zx^ has been included to 
evaluate the 95% confidence interval on the cross sec- 
tions. This procedure is strictly valid only in case of 
purely statistical errors. It is thus not completely sat- 
isfying from this point of view but it should be rather 
conservative however. The gray band on the figure cor- 
responds to this 95% confidence interval obtained for the 
proton production cross section (from the overall proton 
data analysis). 

The obtained confidence intervals on the cross sections 
were then used to compute the uncertainty on the sec- 
ondary particles flux, and ultimately, on the neutrino 
flux. 

For the proton induced kaon production and the pion 
induced pion production the original KMN parametriza- 
tions were used [2^]. These cross sections are less criti- 
cal since their contribution to the final neutrino flux are 
rather small, namely less than 25% and 5% respectively 
up to 30 GeV incident energy. 

E. Particle decay 

In the simulation, the (i/, D, e^ 1 ) spectra from muon de- 
cay, were generated according to the Fermi theory, and 
the muon polarization was taken into account as in the 
previous work by the authors. For the kaon decay, the 
Dalitz plot distribution given in p3j was used. 

F. Variance reduction techniques 

To increase the efficiency of the simulation several tech- 
niques could be used to reduce the variance of the esti- 
mator produced by the simulation j24|. The technique 
used here was the important sampling method which 
consists simply of increasing arbitrarily the incident flux 
of particles, where the latter produce secondaries which 



contribute dominantly to the final differential flux to be 
evaluated, in a particular kinematical domain. 

For instance, in the computation of the high energy 
neutrino flux, the high energy galactic flux of charged 
cosmic rays provides the dominant contribution. The 
natural abundance of this generic flux is extremely small 
however. The above mentioned technique is then used 
to get around this difficulty and the HE CR production 
probability is enhanced by a factor increasing the pro- 
duction yield of the secondaries of interest. The incident 
particle and all the particles produced in this event will 
then simply be weighted with the inverse enhancement 
factor used to boost the particle production. This weight 
is then taken into account when computing the estimator 
from the simulation (see next section). 

Another classical use of this method is for computing 
the flux of a given particle in a given area at a given 
altitude. The primary particle flux generated vertically 
downwards around the zenith of the fiducial virtual de- 
tector area will have a much larger contribution to the 
calculated final flux, compared to CRs produced far away 
from this geographical region. In this case a natural way 
to improve the simulation efficiency is to increase produc- 
tion probability for particles generated near the detection 
region, and to correct for the induced bias by appropri- 
ate statistical weighting of the corresponding events, as 
described above. 

The advantage of the method is that all contributing 
particles are taken into account in the generation process, 
none of them being excluded or neglected. A bad choice 
in the enhancement factor used to increase a certain class 
of primary particles would thus increase the variance of 
the estimators but it would not bias the results. 



G. Particle detection and flux estimation 

Each particle was traced by the program, and his 
history, trajectory parameters and kinematics, were 
recorded. The event file was then analyzed separately 
to generate the various distributions of interest. 

To compare the simulated and the measured flux, a 
virtual surface of detection was defined at the detector 
altitude. It had to be large enough to ensure a large 
counting statistics and a precise estimation of the physics 
observables, but small enough to preserve the accuracy 
of the calculation. In practice the size of area for each 
computation was chosen to provide an accuracy below 
1%. 

The simulated differential particle flux is calculated 
(to -2 • s _1 • sr -1 ■ GeV -1 ) from the sum of the number 
of particles detected in the chosen area times, for each 
particle, the geometrical efficiency in the given particle 
direction of the considered experiment, times the weight 
obtained from the sampling method discussed in section 
III Fl to compensate the enhancement factor. 

This latter number is then divided by the (local) sur- 
face of the particle collection area, the integral of detector 
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acceptance over the full solid angle, energy bin size, and 
equivalent sampling time of the CR flux. 

The equivalent sampling time of CR flux is obtained 
from the total event number generated in the simulation 
run(s) divided by the surface of the generation sphere 
times the integrated (cos0 z weighted) solid angle (tt) 
times the energy integrated flux. 
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FIG. 2: Open squares : Proton flux measured by the BESS 
experiment at mountain altitude 2770 m above sea level, data 
from I25L Light gray band : 95 % confidence interval from the 
simulation. The estimation of the confidence interval includes 
only errors from the production cross sections. The dark gray 
band corresponds to the 95 % confidence interval obtained 
by fitting the cross section on atmospheric data (see section 



III. RESULTS 
A. Proton flux in the atmosphere 

The ability of the simulation to account for the proton 
production in the atmosphere can be probed by compar- 
ing the simulation results with the recent p flux measure- 
ments performed by the CAPRICE experiment between 
sea level and high altitude (~ 40 km) |2u|, and by the 
BESS experiment at mountain altitude [2a ]. 

Figure |21 shows the flux measured by BESS j2f| com- 
pared to the calculated flux. The Bess measurements 
were performed at Mt Norikura, Japan (geographical lo- 
cation 36N,137.5E) at 2770 m of altitude above sea level 
(corresponding to an atmospheric depth of 742 g/cm.2 ) 
in 1999. 

The virtual detector area for this computation was cho- 
sen as the geographical region defined by 0.4 < |Acga/| < 
0.6 rad where Xcgm is the CGM latitude [27]], thus con- 
sisting of two bands (one in each hemisphere) centered 
on plus or minus the experiment CGM latitude (0.509 
rad). 



As in the measurements, the geometrical efficiency is 
set to one for zenithal angular range 9 Z with cos 9 Z > 
0.95. 

Figure |2 shows the confidence interval on the accu- 
racy of the simulated flux (with a 95 % level). The esti- 
mation of the confidence interval includes only the error 
associated to the cross section and were computed us- 
ing the confidence interval associated to the cross section 
parametrization (see section 111 L>|) . The large width of 
the confidence interval in this estimation is due to the 
large number of collisions between the primary particle 
and the detected proton at low altitude. Any uncertainty 
in the proton production cross section is then amplified 
by a large factor (approximatively 8 on the average, see 



by : 

m 



Figure[3]shows the flux measured by CAPRICE for dif- 
ferent altitudes in the atmosphere during the 1998 flight 
from Fort Summer N. M., (34.28 N, 105.14 W, CGM 
latitude 0.76 rad) [2|j and the corresponding simulation 
results. 

The virtual detector area for this computation was cho- 
sen as the region defined by 0.74 < |Acgm| < 0.78 rad. 
The geometrical efficiency used here was taken from [2(| 
to reproduce the detector acceptance. At high altitude, 
the effect of the geomagnetic cutoff is clearly observed in 
all spectra. At this CGM latitude the mean geomagnetic 
cutoff is 6 GeV for protons. Above the cutoff the flux 
is dominated by the primary cosmic component whereas 
below the geomagnetic cutoff the flux consists only of 
secondary protons produced by cosmic rays above the 
cutoff. The primary flux is increasingly absorbed in the 
atmosphere with the decreasing altitude, and the sec- 
ondary component becomes dominant over the entire en- 
ergy range. 

The confidence interval of the calculated flux is seen to 
increase with the decreasing altitude, i.e., with the mean 
number of interactions as it could be expected It 
can be observed that at low altitudes, the uncertainty of 
the simulation results (which includes only the produc- 
tion cross section uncertainties) is much larger than the 
uncertainty from the flux measurements. Therefore, the 
precise proton flux data at low altitudes can be used to 
constrain the parametrization of the secondary particles 
production more tightly and more accurately than the 
available nuclear data. 



B. Muon flux in the atmosphere 

Atmospheric muons are produced in the same de- 
cay chain as neutrinos. Their spectra are thus an es- 
sential ground to probe the reliability of the neutrino 
flux calculated in the same framework. Several exper- 
iments, CAPRICE, HEAT, BESS, have made precise 
flux measurements at various altitudes in the atmosphere 
EHH SIH IHiil which can be used for this purpose. 

Figures 01 and [SI show the calculated positive and neg- 
ative muon flux and the corresponding flux ratios com- 
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FIG. 3: Proton flux measured by the Caprice experiment (open squares) for different altitudes in the atmosphere l2ol compqred 
with the flux estimated from the simulation (solid line). The light gray band shows the 95 % confidence interval from the 
simulation. This estimation of the confidence interval includes only errors from the production cross sections. 
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FIG. 4: Open squares : Muon flux and muon flux ratio measured at Mt Norikura, Japan (location 36 °N, 137.5 °E, altitude 
2770 m above see level corresponding to a atmospheric depth of 742 g/cm2 ). Data from |2Sl| . The light gray band shows 
the 95 % confidence interval from the simulation. The estimation of the confidence interval includes only errors from the 
production cross sections. The dark gray band corresponds to the 95 % confidence interval obtained by fitting the cross section 
on atmospheric data (see text, section PlICI . 



pared with the data measured by the BESS experiment 
at mountain altitude (Mt Norikura, location 36N,137.5E, 
altitude 2770 m) |2^| and at sea level (Tsukuba, location 
36.2 °N,140.1 °E, altitude 30 m) HJ. 

For these two calculations, the virtual detectors were 
chosen as the region defined by 0.4 < \\cgm\ < 0.6 rad 
at the experiment altitude, and with the geometrical effi- 



ciency set to one for zenithal angles 9 Z with cos# z > 0.98. 

The 95 % confidence interval is also shown on the fig- 
ure. It includes the uncertainty originating from the 
parametrization of the secondary protons, neutrons and 
pions production cross sections. Note that for the muon 
charge ratio the uncertainty results only from the uncer- 
tainty on the pion production cross section. 
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FIG. 5: Open squares : Muon Flux and muon flux ratio in Tsukuba, Japan (location 36.2 °N, 140.1 °E, altitude 30 m above 
sea level). The data are from l29t . The light and dark gray bands are defined as in figure |3] 
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FIG. 6: Positive and negative muon flux measured by the Caprice experiment for different altitudes in the atmosphere. The 
solid line shows the calculated flux and the light gray area corresponds to the 95% confidence interval. 
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FIG. 7: Neutrino flux (left) and Neutrino flux ratio (right) at the SK location average over 4ir. The light gray band corresponds 
to the 95 % confidence interval from the simulation. The estimation of the confidence interval includes only errors from the 
production cross sections. The dark gray band represents the 95 % confidence interval obtained by fitting the cross section on 
atmospheric data (see text, section fill CI . 



Figure shows the calculated positive and negative 
muon flux compared with the data measured by the 
CAPRICE experiment during the 1998 flight [3l| at Fort 
Summer. A good agreement is found between data and 
the calculation for the positive and negative fluxes at dif- 
ferent altitude going from the top of atmosphere to the 
ground level. 



C. Use of atmospheric flux measurements to 
constrain hadron production cross sections 

It is clearly seen on the figures [21 to |S] that, at least 
for the low altitude measurements, the statistical preci- 
sion achieved by the atmospheric experiments is signifi- 
cantly better than the precision obtained from the simu- 
lation, taking into account the uncertainty on secondary 
particles production cross section. This results was ob- 
tained although all the relevant available nuclear data 
have been collected to constrain the parametrized cross 
sections used in these calculations. Therefore, the atmo- 
spheric flux measurements are in this sense more accurate 
that the nuclear data. 

In other works on neutrino flux calculations, the calcu- 
lated atmospheric flux of charged particles, used to check 
the reliability of the calculations, were all based on ex- 
isting hadronic Monte Carlo generator 0, [H, |(| , whereas 
in the present approach, in account of the good accuracy 
of the available measurements, the atmospheric data flux 
can be used to constrain the hadronic cross sections to 
ultimately improve the final accuracy of the calculated 
neutrino flux. 

In this work, the proton and muon measurements at 
the Mt Norikura altitude and the muon measurements in 
Tsukuba, from the Bess experiment, were used to con- 



strain both the proton and pion production cross sec- 
tions. The obtained 95 % interval of confidence is shown 
(dark gray) on figure and 151 The reduced \ 2 obtained 
from this fit was 4.7, and the same procedure as described 
in section III Dl was thus used to estimate the confidence 
interval. It must be noted that these fitted cross sec- 
tions may absorb systematic errors originating from other 
sources, like the primary CR fluxes, although the uncer- 
tainty on the latter are much smaller than on the cross 
sections. However, this does not affect the ultimate im- 
provement of the neutrino flux uncertainty, since the de- 
duced interval of confidence accounts for both sources of 
uncertainty. 



D. Neutrino flux at SuperKamiokande Location 

The calculated energy distributions of the atmospheric 
neutrino flux around the SuperKamiokande (SK) detec- 
tor have been computed with the same simulation pro- 
gram. A virtual detector center was defined at the SK 
geographical coordinates (36°N,137°E) with a size of 8° 
in latitude and 18° in longitude (or -900x1600 km 2 ). 
The virtual detector size has been chosen so as not to 
change the calculated flux by more than 1 %. 

Figure shows the average flux over the full 4ir range 
of solid angle for + and v e + D e (left) and the flux 
ratio i |r and |^ (right). The light gray area rep- 

resents the 95 % confidence level corresponding to the 
uncertainty due to the particle production cross sections 
inaccuracies. The uncertainty for the absolute neutrino 
flux is of the order of 10 % and, as expected, these 
error contributions are largely reduced for the ratio of 
flux I s (insensitive to proton/neutron production cross 
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FIG. 8: Zenith angle distributions and flavor ratio (y^ + P M )/(i/ e + i> e ) of the calculated neutrino flux around the Super- 
Kamiokande detector for five energy bin from 0.1 to 30 Gev. The light and dark gray bands are defined as in Fig. (see text, 
section Ml Cll . 



sections) and for ^- (insensitive to proton/neutron pro- 
duction cross sections but also to pion production cross 
sections for low energy neutrinos) and vanished for the 
^''^p'' flux ratio (insensitive to proton/neutron and pion 
production cross sections). The dark gray band corre- 
spond to the 95 % confidence interval obtained by fitting 
the cross section on atmospheric data discussed in the 
section till CI Using this method the uncertainty for the 
absolute flux is reduced to the order of 3 %. 

The zenith angle distribution of the flux and of the 
flavor ratios are shown on figure |S1 left for 5 energy bins 
between 0.1 and 30 GcV. At low energy the zenith an- 
gle distribution displays an enhancement for directions 
close to horizontal, in qualitative agreement with |34j . 
This enhancement originates from the isotropisation of 
the neutrino flux resulting from cosmic ray showers. It is 
expected to be maximum for a purely isotropic neutrino 
emission and nil for a perfectly collimatcd neutrino pro- 
duction. The maximum is then naturally obtained for the 



lowest energy bin (about twice the downward or upward 
going flux out of the horizontal plane). It is decreas- 
ing with the increasing energy. At low energy the zenith 
angle distribution for downward and upward neutrinos, 
outside the cos = peak, is approximatively flat, this 
isotropy naturally originating from the CR isotropy. The 
small structures observed arise from geomagnetic effects 
on the primary flux. At high energies the distributions 
become sensitive to the muon decay probability which 
depends on the muon pathlcngth in the atmosphere and 
on its energy. As expected the | cos 9 Z | dependence is be- 
coming steeper with the increasing energy. The effect is 
enhanced for electronic neutrino which originate exclu- 
sively from muon decay. This is shown on figure El right 
where the flavor ratio (f^ + P^) / \v e + D e ) is seen to depend 
on the zenith angle, especially for high energy neutrinos, 
for the same reason as given previously. 

On these figures, the light gray area represents the 95% 
confidence level corresponding to the particle production 
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FIG. 9: Azimuth angle distribution of the neutrino flux (left) and flavor ratios, v e /v e and v^fv^ (right) at the SK detector 
location for five energy bins from 0.1 to 30 Gev. The light and dark gray bands are defined as in Fig. (see text, section lillCft . 



uncertainty. Here again we see that the uncertainty is 
reduced at high energy, this is mainly due to the de- 
crease of the mean number of interactions between the 
primary cosmic ray impinging the atmosphere and the fi- 
nal neutrino detection. As in figured f°r the flavor ratio 
(ia, + 9n)/{v e + v e ) the error contributions from particle 
production uncertainty vanished. The dark gray band 
represents the 95 % confidence interval obtained by fit- 
ting the cross section on atmospheric data as described 
in section UlI CI 

The azimuth angle distribution of the flux and of the 
flavor ratios are shown on figure for 5 energy bins be- 
tween 0.1 and 30 GcV. On this figure northward and 
westward going particles correspond to azimuth angles 
0, | respectively. The non-flatness of the low energy bin 
originates from the geomagnetic cutoff on primary cosmic 
rays [3{| . The east- west (EW) effect of the primary flux 
(more eastward going than westward going cosmic rays) 
results in a EW asymmetry in the produced secondary 
particles and then in the neutrino flux at low energy. 



With the increasing energy, neutrinos are increasingly 
produced in atmospheric showers initiated by higher en- 
ergy cosmic rays which are less and less sensitive to the 
geomagnetic cutoff (the cosmic ray flux becomes isotropic 
at a rigidity of 60 GV). This results in the simulated 
azimuthal distribution becoming flat for the last energy 
bin in figure |5Jlcft (neutrino energy between 10 Gev and 
30 Gev). 

For the azimuth angle distributions of the flavor ratios 
(figure EJright), it can be noted that the v e /v e ratio is 
strongly EW asymmetric while the opposite is observed 
for the v^lvp ratio. The (y^ + i7^)/{v e + i7 e ) ratio (not 
shown) is found to be almost structureless at all latitudes. 
This additional EW asymmetry featured by the u e /P e 
and v^jv^ ratio can be explained by the muon bending 
in the geomagnetic field [35j. In the production chain 
p — ► 7r + — > /-t + f M — > VeVfi, the muon will propagate with 
the same bending as the proton one. An enhancement 
of EW effect originating from the geomagnetic cutoff is 
then expected for the neutrino v e and produced in the 
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(i + decay. In the p — > ir~ — > pTv^ — > VyV e production 
chain, the muon will propagate with the opposite bend- 
ing and a reduction of the EW asymmetry is expected 
for the \i~ decay products and v e . This explains the 
EW features seen in the v e /D e and the Vfi/v^ ratio. In 
addition, because muonic neutrinos are also produced in 
the pion decay, the EW asymmetry is found to be higher 
for v e /v e than for v u .lv n - These results confirm the qual- 
itative prediction from |35j : 

Ap e < A Vfi < Afj^ < A Ve 

where A stands for the neutrino EW asymmetry. 

In figure [U] the light and dark gray area represents 
respectively the 95 % confidence level corresponding to 
the particle production cross section uncertainty and the 
95 % confidence interval obtained by fitting the cross sec- 
tion on atmospheric data as described in section IIII CI 
The uncertainty are largely reduced for the ratio of flux 
|^ (insensitive to the proton/neutron production cross 

sections) and for ^ (insensitive to proton/neutron pro- 
duction cross sections but also to pion production cross 
sections for low energy). 

IV. SUMMARY AND CONCLUSION 

A full three-dimensional simulation has been used to 
compute the secondary particles production in the atmo- 
sphere. In these calculations the statistics of the simu- 
lated sample has been significantly increase with respect 



to our previous work, resulting in a negligible contribu- 
tion of the statistical uncertainty. The dominant error 
then appeared to come from the uncertainty on the par- 
ticle production cross sections. In this work this system- 
atic uncertainty has been studied quantitatively. 

The simulation has been used to reproduce the atmo- 
spheric measurements of the proton and muon flux, in 
order to test the reliability of the calculations. These 
atmospheric measurements, in account of their good ac- 
curacy, were also used to constrain the secondary parti- 
cle production cross sections, which experimental values 
from direct measurements on accelerators were of lesser 
accuracy. 

The results for the absolute neutrino fluxes and 
their zenithal and azimuthal angular distributions to- 
gether with the flavor ratio of these flux at the Supcr- 
Kamiokande location, have been presented. The high 
statistics accumulated allowed to investigate the detailed 
structure of these distribution and to discuss them. The 
uncertainty of the calculations has been estimated. It 
has been shown first, that this uncertainty is at the level 
of 10 % for the absolute flux, and second that it can be 
reduced to the level of 3 % by using the atmospheric data 
on the muon and proton flux to constrain the production 
cross sections. These results provide an accurate basis 
for the comprarison to future experimental data. 

Acknowledgments 

The author is very grateful to Michel Buenerd without 
whom this work would not have been initiates and for his 
precious help in writing the manuscript. 



[1] Super-Kamiokande Collaboration, Y. Fukuda, et al., 
Phys. Rev. Lett. 81, 1562(1998), B433, 9(1998), B436, 
33(1998). 

[2] The Soudan 2 Collaboration, W. W. M. Allison, et al., 
Phys. Lett. B449, 137(1999). W. W. M. Allison, et al., 
Phys. Lett. B391, 491(1997). 

[3] J. Wentz, Phys. Rev. D 67, 073020 (2003). 

[4] M. Honda et al., Phys. Rev. D 70, 043008(2004). 

[5] G. Battistoni et al. Astropart. Phys. 19, 269(2003). 

[6] CD. Barr et al., Phys. Rev. D 70, 023006(2004). 

[7] L. Derome et al., Phys. Lett. B489, 1(2000). 

[8] L. Derome, M. Buenerd, and Yong Liu, Phys. Lett. B515, 
1(2001). 

[9] L. Derome, and M. Buenerd, Phys. Lett. B521, 
139(2001). 

[10] Y. Liu, L. Derome, M. Buenerd, Phys. Rev. D 67 
073022(2003). 

[11] B. Baret, L. Derome, and M. Buenerd, Phys. Rev. D68, 
053009(2003). 

[12] Y. Liu, L. Derome and M. Buenerd, 27th Int. Cosm. Ray 

Conf., Hamburg, Aug 8-13, 2001. 
[13] The AMS collaboration, J. Alcaraz et al., Phys. Lett. 

B472, 215(2000), Phys. Lett. B490, 27(2000), Phys. Lett. 

B494, 19(2000). 
[14] B. Wiebel-Sooth et al., Astronomy and Astrophysics 330, 



389(1998) 

[15] J.S. Perko, Astro. Astrophys. 184(1984)119; see M.S. 
Potgieter, Proc. ICRC Calgary, 1993, p213, for a recent 
review of the subject. 

[16] J. Engel, T. K. Gaisser, T. Stanev and P. Lipari, Phys. 
Rev. D46, 5013 (1992). 

[17] J.V. Allaby et al., CERN Yellow report 70-12, (1970). 

[18] A. E. Hedin, J. Geophys. Res. 96, 1159(1991). 

MSISE model available at 

http : //modelweb .gsf c .nasa.gov/ atmos/ 

[19] the IGRF models are available at 
http : //www .ngdc .noaa.gov/IAGA/vmod/home . html 

[20] W. H. Press, S. A. Teukolsky, W. T. Vetterling, Numer- 
ical Recipes in C++: The Art of Scientific Computing, 
Cambridge University Press (2002) 

[21] L. Derome and M. Buenerd, Proc. 29th ICRC, Pune (In- 
dia), 2005, HE.2.1 

[22] A. N. Kalinovsky, N. V. Mokhov and Y. P. Nikitin, Pas- 
sage Of High-Energy Particles Through Matter, NEW 
YORK, USA: AIP 1989. 

[23] S. Eidelman et al. (Particle Data Group), Physics Letters 
B592, 1(2004). 

[24] A. Haghighat and JC. Wagner, Prog. Nucl. Energy 42, 
25(2003). 

[25] T. Sanuki et al., Phys. Lett. B577, 10(2003). 



11 



[26] E. Mocchiutti, PhD thesis, Royal Institute of Technology, 

Stockholm, 2003 
[27] Gustafsson, G., N. E. Papitashvili, and V. O. Papi- 

tashvili, J. Atmos. Terr. Phys., 54, 1609(1992). 
[28] BESS Collaboration, T. Sanuki et al., Phys. Lett. B541, 

234(2002). 

[29] BESS Collaboration, S. Haino et al. Phys. Lett. B594, 
35(2004). 

[30] BESS Collaboration, M. Motoki et al., Astropart. Phys. 
19, 113(2003). 



[31] CAPRICE Collaboration, M. Boezio et al., Phys. Rev. D 

67,072003 (2003). 
[32] CAPRICE Collaboration, J. Kremer, et al., Phys. Rev. 

Lett. 83, 4241(1999). 
[33] HEAT Collaboration, S. Coutu, et al, Phys. Rev. D62, 

032001(2000). 
[34] P. Lipari, Astropart. Phys. 14, 153(2000) 
[35] P. Lipari, Astropart. Phys. 14, 171(2000) 



